#!/Pub/Apps/bin/Rscript
suppressPackageStartupMessages(library(argparse))

wrp <- "转化生成loom文件，需要两个参数，第一个为seurat对象路径(rds or rdata)，第二个为od输出路径。可选conuts或者data。"
parser <- ArgumentParser(
    prog='seurat2loom_exprMat.r',
    description = wrp,
    formatter_class = "argparse.RawTextHelpFormatter"
    # usage = "seurat2loom_exprMat.r -h; seurat2loom_exprMat.r seurat_obj od"
)

parser$add_argument("seurat_obj", nargs = 1, help = "seurat对象路径。[rds or rdata]")
parser$add_argument("od", nargs = 1, help = "结果输出路径")
parser$add_argument("-t", "--type_of_data",
    action = "store",
    default = "data", metavar = "",
    help = "分析所用数据类型，貌似counts或者data都可以。默认data"
)

args_in_process <- parser$parse_args()


if (file.exists(stringr::str_glue("{args_in_process$od}/seurat_obj_expression.loom"))) {
    cli::cli_alert_warning("\nseurat_obj_expression.loom已生成，exprMat转化lomm跳过\n")
    q()
}

dir.create(args_in_process$od,recursive = T,showWarnings = F)
Sys.setenv(HDF5_USE_FILE_LOCKING = "FALSE")

suppressMessages(library(Seurat))
suppressMessages(library(SCopeLoomR))
suppressMessages(library(stringr))

cli::cli_inform(">> 1.1 Load Seurat obj [{Sys.time()}]\n")

fn <- basename(args_in_process$seurat_obj)
if (grepl("rds", fn, ignore.case = T)) {
    ss_tmp <- readRDS(args_in_process$seurat_obj)
}
if (grepl("rdata", fn, ignore.case = T)) {
    tmp <- load(args_in_process$seurat_obj)
    ss_tmp <- get0(tmp)
}

f <- function(){
    cli::cli_alert_info("Seurat info:\n")
    print(ss_tmp)
    cat('\n')
}

# http://htmlpreview.github.io/?https://github.com/aertslab/SCENIC/blob/master/inst/doc/SCENIC_Setup.html
# Overall, SCENIC is quite robust to this choice, we have applied SCENIC to datasets using raw (logged) UMI counts, normalized UMI counts, and TPM and they all provided reliable results (see Aibar et al. (2017)).

exprMat = switch(args_in_process$type_of_data,
    data = {
        if (startsWith(as.character(utils::packageVersion("Seurat")), "5")) {
            exprMat <- ss_tmp@assays$RNA$data
        } else {
            exprMat <- ss_tmp@assays$RNA@data
        }
    },
    counts = {
        if (startsWith(as.character(utils::packageVersion("Seurat")), "5")) {
            exprMat <- ss_tmp@assays$RNA$counts
        } else {
            exprMat <- ss_tmp@assays$RNA@counts
        }
    }
)

cellInfo <- ss_tmp@meta.data

loci1 <- which(rowSums(exprMat) > 1 * .01 * ncol(exprMat))
exprMat_filter <- exprMat[loci1, ]

# function from https://bookdown.org/ytliu13207/SingleCellMultiOmicsDataAnalysis/scenic.html#retrieve-auc-scores-per-cell
# 貌似不需要添加cell信息到
add_cell_annotation <- function(loom, cellAnnotation) {
    cellAnnotation <- data.frame(cellAnnotation)
    if (any(c("nGene", "nUMI") %in% colnames(cellAnnotation))) {
        warning("Columns 'nGene' and 'nUMI' will not be added as annotations to the loom file.")
        cellAnnotation <- cellAnnotation[, colnames(cellAnnotation) != "nGene", drop = FALSE]
        cellAnnotation <- cellAnnotation[, colnames(cellAnnotation) != "nUMI", drop = FALSE]
    }

    if (ncol(cellAnnotation) <= 0) stop("The cell annotation contains no columns")
    if (!all(get_cell_ids(loom) %in% rownames(cellAnnotation))) stop("Cell IDs are missing in the annotation")

    cellAnnotation <- cellAnnotation[get_cell_ids(loom), , drop = FALSE]
    # Add annotation
    for (cn in colnames(cellAnnotation))
    {
        add_col_attr(loom = loom, key = cn, value = cellAnnotation[, cn])
    }

    invisible(loom)
}

cli::cli_inform(">> 1.2 Converting loom [{Sys.time()}]\n")
loom <- SCopeLoomR::build_loom(
    stringr::str_glue("{args_in_process$od}/seurat_obj_expression.loom"),
    dgem = exprMat_filter
)
loom <- add_cell_annotation(loom, cellInfo)
close_loom(loom)


if (FALSE) {
    # 结果分析示例
    library(SCopeLoomR)
    library(SCENIC)

    loom <- SCopeLoomR::open_loom("out/4.sc_houxu/十.3/diercichangshi/out_SCENIC.loom")
    regulons_incidMat <- get_regulons(loom, "Regulons")

    regulons_incidMat <- get_regulons(loom, column.attr.name = "Regulons")

    regulons_incidMat[1:4, 1:4]
    regulons <- regulonsToGeneLists(regulons_incidMat)

    regulonAUC <- get_regulons_AUC(loom, column.attr.name = "RegulonsAUC")
    regulonAucThresholds <- get_regulon_thresholds(loom)
    tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))])

    embeddings <- get_embeddings(loom)
    close_loom(loom)

    saveRDS(list("regulons" = regulons, "regulonAUC" = regulonAUC),
        file = "out/4.sc_houxu/十.3/diercichangshi/ss_t_FilterLQCells-pyscenic-output-AUC.rds"
    )

    # 一些可视化参考 http://www.360doc.com/content/22/1104/09/74726315_1054515292.shtml
}

